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Understanding how laser light scatters from realistic mirror surfaces is crucial for the design, com¬ 
missioning and operation of precision interferometers, such as the current and next generation of 
gravitational-wave detectors. Numerical simulations are indispensable tools for this task but their 
utility can in practice be limited by the computational cost of describing the scattering process. 
In this paper we present an efficient method to signihcantly reduce the computational cost of op¬ 
tical simulations that incorporate scattering. This is accomplished by constructing a near optimal 
representation of the complex, multi-parameter 2D overlap integrals that describe the scattering 
process (referred to as a redueed order quadrature). We demonstrate our technique by simulating a 
near-unstable Fabry-Perot cavity and its control signals using similar optics to those installed in one 
of the LIGO gravitational-wave detectors. We show that using reduced order quadrature reduces 
the computational time of the numerical simulation from days to minutes (a speed-up of 2750 x) 
whilst incurring negligible errors. This signihcantly increases the feasibility of modelling interfer¬ 
ometers with realistic imperfections to overcome current limits in state-of-the-art optical systems. 
Whilst we focus on the Hermite-Gaussian basis for describing the scattering of the optical helds, our 
method is generic and could be applied with any suitable basis. An implementation of this reduced 
order quadrature method is provided in the open source interferometer simulation software Finesse. 


I. INTRODUCTION 

Laser interferometers have long been an exceptional 
tool for enabling high precision measurements. With 
ever increasing demands on their performance, new tech¬ 
niques and tools have been developed to design and build 
the next-generation of instruments. This has especially 
been true in the development of gravitational-wave de¬ 
tectors over the last several decades [1-4]. Such ground- 
based gravitational-wave detectors are based on a Michel- 
son interferometer and are enhanced with Fabry-Perot 
cavities. Detecting gravitational waves is still one of the 
major challenges in experimental physics, and the inter¬ 
ferometers used include numerous new optical technolo¬ 
gies to reach unprecedented displacement sensitivities be¬ 
yond 10“^^ m/v/Sz. 

Some of these detectors are currently being upgraded 
to have a ten-fold increase in sensitivity using a much 
higher circulating power [5, 6]. To achieve their target 
performance the detectors undergo several years of com¬ 
missioning, during which the interferometers are carefully 
tested and improved towards their designed operational 
state. Numerical simulations are important tools to di¬ 
agnose causes of any unexpected behaviour seen during 
commissioning; to suggest solutions to potential prob¬ 
lems and for advising the design of detector upgrades. 
Hence, there is a long history of developing and using 
dedicated optical simulation tools for the commissioning 
and design of gravitational-wave detectors [7-10]. 

One of the key aspects for the current instruments is 
the high circulating laser power, up to hundreds of kilo- 
Watts, required for a broadband reduction of shot-noise. 
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It has been recognised for some time that the thermal 
deformations of the optics due to spurious absorption 
can degrade the performance of the interferometers [11]. 
Numerical models have been used extensively in the in¬ 
vestigation of such problems and in the development of 
mitigating solutions (for example [12, 13]). Thermally 
induced distortions and other effects related to the laser 
beam shape are still limiting factors of the instruments 
today and are concerns for the design of future detec¬ 
tors [14]. Furthermore, similar effects can limit the per¬ 
formance of other optical precision measurements such 
as optical clocks [15] or the optical readout of atomic 
systems [16]. Mitigation strategies for beam shape dis¬ 
tortions in complex interferometers are actively being de¬ 
veloped and require accurate numerical models for their 
design and development. 

Initially, the simulation tools for investigating dis¬ 
torted beams used a grid-based field description. Beam 
distortions can also be modeled effectively using an ex¬ 
pansion into spatial cavity eigenmodes [17], such as 
Hermite-Gauss modes. The interaction of the beam 
shape with a distorted optical surface often requires the 
computation of a scattering matrix based on measured or 
simulated profiles of the distorted surface. This is always 
true for mode-based simulations programs but is also re¬ 
quired for grid-based codes when specific shapes of the 
beam are important, for example, for the investigations 
of parametric instabilities [? ? ]. If this matrix has to be 
re-generated, for example when the effects of a change of 
a surface shape is being investigated, this element of the 
computation can dominate the total time required for the 
entire simulation. A prominent example is that when the 
circulating laser power within the LIGO interferometers 
thermally warps the mirror surfaces changing the shape 
of the laser beams and requiring a re-calculation of many 
scattering matrices. Including this effect can increase the 
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computation time from minutes to days. 

Some of us are providing numerical simulation support 
for the commissioning of the LIGO interferometers [18]. 
We use our own simulation tool Finesse [19] and are 
maintaining parameter files for the detectors [20]. Com¬ 
missioning tackles the unexpected behaviour of the in¬ 
terferometers and must take into account the sometimes 
rapid progress of the experimental setup. Therefore sup¬ 
port provided with numerical models must fulfil two cri¬ 
teria: a) we must be able to accurately model the current 
experimental setup in the presence of distortions and de¬ 
viations from the design and b) we must be able to pro¬ 
vide a quick response to new questions to inform the 
management of the activities on site in real time. Fi¬ 
nesse is a frequency-domain tool, using Hermite-Gauss 
modes to describe beam distortions and is thus ideally 
suited as a rapid and accurate tool. 

Our investigations with numerical tools typically con¬ 
sist of a sequence of different subtasks, sometimes using 
different tools, alternating with an expert review of in¬ 
termediate or preliminary results. This is a very different 
pattern of tasks to those that benefit from a computer 
cluster or super computer. Instead, our work requires 
lightweight and flexible tools with computing times up 
to minutes or hours. Because of this, strategies to ame¬ 
liorate the run time of simulations are of high importance 
to provide fast diagnosis of unexpected behaviour; to al¬ 
low the parameter space of the simulations to be probed 
exhaustively; to improve the resolution of simulations at 
a fixed run-time, and to allow simulations to be run on 
less powerful and cheaper hardware. 

In this paper we present a new approach that reduces 
the computational time of simulations based on modal 
models by several orders of magnitude. We specifically 
target the computational cost of computing scattering 
matrices for optical simulations. Our approach is based 
on a near-optimal formulation of the integrals required 
to compute the scattering matrices, known as a reduced 
order quadrature [21] (ROQ). The reduced order quadra¬ 
ture has already been applied in the context of astro¬ 
nomical data analysis with LIGO [22] where the repeated 
computation of quantities similar to the scattering matrix 
dominate the run time of the analysis codes. Crucially, 
the reduced order quadrature is designed to provide huge 
improvements to computational efficiency whilst main¬ 
taining computational precision. 

The ROQ can be regarded as a type of near-optimal, 
application specific, downsampling of the integrands 
needed to compute the integrals for the scattering ma¬ 
trices [21]. It is analogous to Gaussian quadrature, but 
whereas Gaussian quadrature is designed to provide ex¬ 
act results for polynomials of a certain degree, the ROQ 
produces nearly-exact results for arbitrary parametric 
functions. Importantly, we are able to place tight er¬ 
ror bounds on the accuracy of the ROQ for a particular 
application [21] making it an ideal technique to speed up 
costly integrals. It exploits an offline/online methodol¬ 
ogy in which we recast the expensive integrals used to 


compute scattering matrices into a more computation¬ 
ally efficient form in the “offline” stage. This is then 
used for the rapid “online” evaluation of the scattering 
matrices. The offline stage can itself be computationally 
expensive, however it need only be performed once and 
is easily parallelised. The data computed in the offline 
stage—that is needed by the ROQ—can be stored and 
shared for particular use cases in the online stage so that 
the offline cost does not need to be factored in at run 
time. 

We derive the algorithm in a general form and report 
on the implementation and performance of this method 
in an example task for the LIGO interferometers. The 
implementation of the method described in this article is 
available as open source as part of the Finesse source 
code and the Python based package Pykat [? ], which 
will also contain the offline computed data to enable oth¬ 
ers to model Advanced LIGO like arm cavities. Our par¬ 
ticular implementation here is used to provide a simple, 
real-word example. However, the algorithm can be eas¬ 
ily implemented in other types of simulation tools, for 
example, time domain simulations or grid based tools 
(also known as FFT simulations) that compare beam 
shapes. In all cases our algorithm can significantly reduce 
the computation time for evaluating overlap integrals of 
Gaussian modes with numerical data. 

The paper is outlined as follows: In section II we give 
an overview of the paraxial description of the optical 
eigenmodes and scattering into higher order modes. In 
Section III we provide the mathematical background and 
algorithm for producing the ROQ. Section III heavily re¬ 
lies on an additional mathematical technique known as 
the “empirical interpolation method” [23]. We assume 
no prior knowledge of this and provide the main details 
and results necessary for the ROQ. Section IV then high¬ 
lights an exemplary case to demonstrate our method for 
modelling near-unstable optical cavities. Finally in sec¬ 
tion V the computational performance of our method is 
analysed. 


II. HIGHER-ORDER OPTICAL MODES 

Gravitational wave detectors are constructed of multi¬ 
ple optical cavities based on a Michelson interferometer. 
The circulating laser beams in such an optical setup is 
well described by the the paraxial Gaussian eigenmodes 
of a spherical cavity; an efficient basis for describing the 
spatial properties of a laser beam in the transverse plane 
to the propagation axis [24]. The fundamental Gaussian 
mode is described in cartesian coordinates by: 
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Wx and Wy are the beam spot sizes in the x and y di¬ 
rections, k is the wavenumber of the laser light and 
q = {qx^Qy} are the complex beam parameters in the 
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X and y directions. The shape of the Gaussian mode is 
fully defined by the wavelength of the light A and the 
beam parameter: 

Qx — ^^R,x — T 1 ^ (2) 

where 2 : 3 ^ is the distance from the waist, zr^x is the 
Rayleigh range, is the size of the waist and A = 

1064nm is the wavelength of the NdiYAG laser used in 
current GW detectors. The same set of parameters exists 

for Qy. 

Any perturbations in the beam’s spatial profile from 
this fundamental Gaussian can be described by the ad¬ 
dition of higher-order Gaussian modes (HOMs). In this 
paper we discuss in particular the cartesian orthogonal 
basis of Hermite-Gaussian (HG) modes [24], however our 
method is applicable for any other suitable basis. The 
complex transverse spatial amplitude of these HG modes 
is given by: 


z{x^y). Thus, both the amplitude and phase of the beam 
will be affected and A(x, y) = a{x^ y)e^‘^^ z{x,y) ^ exam¬ 

ple of the measured surface height variations present on 
LIGO test mass mirrors can be seen in figure 1 [25]. The 
mode content of the outgoing beam E{x^y; q) is com¬ 
puted by projecting E' into the outgoing beam basis q. 
For any incoming HOM Un>m' fhe amount of outgoing 
Unm can be computed via an overlap integral, this com¬ 
plex value is known as a coupling coefficient: 

^nm,n'mfQx-)Qx')Qy")Qy ''>^) ~ 

/ / IC{Xx] x)A{x, y)IC{Xy] y) dx dy , (5) 

J J —00 

where the integral kernels ]C{Xx]x) and IC{Xy;y) are 
given by 

K:{\oc\x) =u*^(x,q^)un'{x,q'^), ( 6 ) 

=u*m{y,qy)um'{y,qy), (7) 
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where n defines the order of the Hermite polynomials 
in the x axis and m for the y. The order of the optical 
mode is O = n + m and individual modes are typically 
referred to as TEM^m. A laser field with a single optical 
frequency component uj can be expanded into a beam 
basis whose shape is described by q as: 


n-\-m<0 

max 

E{x,y,t;q)= 2 /; q)e‘“‘ (4) 

n=0,m=0 


where a^m is a complex value describing the amplitude 
and phase of a mode TEM^m and Omax is the maximum 
order of modes included in the expansion. 


A. Scattering into higher-order modes 

When a field interacts with an optical component its 
mode content is typically changed. Here we define scat¬ 
tering as the relationship between the mode content of 
the outgoing beams d with a beam shape q, and the 
mode content of the incoming beam d' described with 
q'. Mathematically this is simply d = kd' where k is 
known as the scattering matrix. Now consider the spa¬ 
tial profile of a beam reflected from on an imperfect optic 
E'{x^y]q') = A(x,^)Ein(x,q'), where E[n is the inci¬ 
dent beam and A{x,y) is complex function describing the 
perturbation it has undergone. Eor example, on reflec¬ 
tion a beam will be clipped by the finite size of the mirror 
a{x^y) and reflected from a surface with height variations 


and the parameter vectors are given by = 

{n^n' ^Qx^q'x) and Xy = (m, m', g^). There are two 

general cases when computing (5): g 7 ^ g' which we refer 
to as mode-mismatched and g = g' as mode-matched. 

Gomputing the scattering matrix k requires evaluat¬ 
ing the integral (5) for each of its elements. If couplings 
between modes up to and including order O are con¬ 
sidered then the number of elements in k is Nj^{0) = 
{O^ + 60^ - 1 - 130^ - 1 - 120 + 4)/4 and the computational 
cost of evaluating this many integrals can be very expen¬ 
sive. In our experience [18, 26] a typical LIGO simulation 
task involving HOMs can be performed with O = 6 — 10 
while some cases, such as those that include strong ther¬ 
mal distortions or clipping, a higher maximum order is 
required. 

In simple cases where A{x^y) = 1 or A{x^y) repre¬ 
sents a tilted surface, analytical results are available for 
both mode matched and mismatched cases [27, 28]. In 
general however A{x,y) is of no particular form and the 
integral in (5) must be evaluated numerically. It is pos¬ 
sible to split multiple distortions into separate scattering 
matrices A{x,y) ^ A{x,y) B{x,y) and the coupling co¬ 
efficients become a product of two separate matrices: 

^nm,n'm' (Q? Q 5 ^ -^) — 

00 

^nm,nm(Q5 Q 5 ^fifh,n'm' (4q';-B) ( 8 ) 

n,m=0 

where q is an expansion beam parameter which we are 
free to choose. Thus our scattering matrix becomes 
^(q>q') = ^A(q,q)^B(q,q')- By choosing q = q or q' 
we can set the mode-mismatching to be in either one ma¬ 
trix or the other. This is ideal as a mode-matched k is 
a Hermitian matrix whose symmetry can be exploited to 
only compute one half of the matrix. By ensuring that 
this matrix also contains any distortions that require nu¬ 
merical integration the computational cost can be nearly 
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Figure 1: Measured surface distortions for the mirrors currently installed in the Livingston LIGO site (shown here are the 
distortions of the test masses before they were coated). ETM09 and ITM08 are installed in the x-arm and ETM07 and 
ITM04 in the y-arm [25]. Theses measurements have been processed to remove the overall mirror curvatures, offset and tilts. 


halved. It is then possible to benefit from the fast an¬ 
alytic solutions to (5) to account for mode-mismatching 
in the other matrix. 


III. EFFICIENTLY COMPUTING SCATTERING 
MATRICES: INTEGRATION BY 
INTERPOLATION 

For a discretely sampled mirror map with L sample 
points in both the x and y directions, the coupling coeffi¬ 
cient (5) can be approximated using a composite Newton- 
Cotes quadrature rule: 

^nm,n'm' {<lx,qx,qy,qy\A) 

L L 

-EE yi) (9) 

k=l1=1 

where LE is an L x L matrix describing the 2D composite 
Newton-Cotes quadrature weights over the area of the 
map. The matrix is found by taking the outer product of 
the ID composite Newton-Cotes quadrature weights [29] 
in both X and y directions. There are terms in the 
double sum (9). When is large, as in the cases of 
interest for this paper, there are two major bottlenecks: 

(i) evaluation of the kernel at each discrete Xk^yk and, 

(ii) evaluation of the double sum. With a set of M ^ L 
basis elements that accurately spans the kernel space, it 
is possible to replace the double sum (9) with a reduced 
order quadrature (ROQ) rule (19) containing only 
terms, reducing the overall cost of the by a factor of ^ 

provided the kernel can be directly evaluated. 

The reduced order quadrature scheme is implemented 
in three steps. The first two are carried out offline, 
while the final, mirror-map-dependent step is performed 
in preparation for the simulations; once per map. The 
steps are as follows: Step 1 - Construct a reduced basis 
(offline); a set of M basis elements whose span describes 


the kernel space. Step 2 - Construct an interpolant us¬ 
ing the basis (offline) by requiring it to exactly match any 
kernel at M carefully chosen spatial subsamples 
[30] (and similarly for y). Step 3 - Use the interpolant 
to replace the inner product evaluations in (9) with the 
ROQ (19) (online). 


A. The Empirical Interpolation Method 

The empirical interpolation method is an efficient tech¬ 
nique performing this offline/online procedure and has 
been demonstrated in the context of astronomical data 
analysis with LIGO [22]. Provided the kernels vary 
smoothly with over x and Xy over y then there exists a 
set of kernels at judiciously chosen parameter values that 
represent any kernel - and hence any integral (5) - for an 
arbitrary parameter value. This set of kernels constitutes 
the reduced basis: Given any parameter value or Xy 
we can And the best approximation to the kernel at A^, 
or Xy as linear combination of the reduced basis. 

The ability to exploit the reduced basis to quickly eval¬ 
uate (5) depends on being able to And an affine param¬ 
eterization of the integral kernels. In general, the ker¬ 
nels do not admit such a parameterization. However, 
the empirical interpolation method finds a near-optimal 
affine approximation whose accuracy is bounded by the 
accuracy of the reduced basis [21]. This affine approx¬ 
imation is called the empirical interpolant. The spatial 
integrals over dxdy in (5) will only depend on the re¬ 
duced basis (and hence only have to be computed once 
for a given mirror map) and the parameter variation is 
handled by the empirical interpolant at a reduced com¬ 
putational cost. 

The empirical interpolation method exploits the of¬ 
fline/online computational concept where we decompose 
the problem into a (possibly very) expensive offline part 
which affords a cheap online part. In this case, the ex¬ 
pensive offline part is in finding the reduced basis and 
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constructing the empirical interpolant. Once the empiri¬ 
cal interpolant is found then we use it for the fast online 
evaluation of (5). One of the main reasons why the em¬ 
pirical interpolant is used for fast online evaluation of 
(5) is due to its desirable error properties that makes it 
superior to other interpolation methods, such as polyno¬ 
mial interpolation. In addition, the empirical interpolant 
avoids many of the pitfalls of high-dimensional interpo¬ 
lation that we would otherwise encounter (see, e.g. [31]). 

B. Affine Parameterization 

We would like the kernel to be separable in the mode 
parameters (A^,, Xy) and spatial position (x^y). For these 
reasons we will look for a representation of the kernel that 
has the following form: 

K{Xx-,x) = a{Xx)f{x), 

= a{Xy) f{y). (10) 

The functions a and / are the same irrespective of 
whether the kernel is a function of x or ^ due to the sym¬ 
metries of the Hermite Gauss modes. Using the affine 
parameterization, the coupling coefficient (5) is: 

^nrn,n'rn'{Q x") Qx") ~ ^ 

[f f*ix)A{x,y)f{y)dxdy, (11) 

This affine parameterization thus allows us to compute all 
the parameter-dependent pieces efficiently in the online 
procedure as all the x — y integrals are performed only 
once for a given mirror map. In general the kernel will 
not admit an exact affine decomposition as in (10). Using 
the EIM, the approximation to the kernels will have the 
form: 

fC{X::,;x) Ci{Xx)ei{x), (12) 

IC{Xy-, y)^Yl ■ 

i 

The sum is over the reduced basis elements and coef¬ 
ficients Ci that contain the parameter dependence. 

Given a basis ei(x), the Ci{Xx) in (12) are the solu¬ 
tions to the M-point interpolation problem whereby we 
require the interpolant to be exactly equal to the kernel 
at any parameter value A^, at a set of interpolation nodes 

M M 

lC{X,-,Xj) = J2ci{XMXj) =J2^jiCi{X,), (13) 

i=l i=l 

where the matrix V is given by 



( ei(W) 

eHXi) ■■ 

• e^(Xi) \ 




eHX2) ■■ 

• e^(X2) 


u = 

e\X^) 

e^Xs) ■■ 

■ e^{Xs) 

(14) 


\e\XM) 

cHXm) ■■ 

■ e^{XM) ) 



Thus we have: 

M 

Ci{X,) = Y,{V~")i^l^{K\Xj). (15) 

i=i 

Substituting (15) into (12), the empirical interpolant 
is: 

M 

Im[JC]{X^;x) =Y,lC{X,-,Xj)Bj{x) (16) 

i=i 

where: 

M 

By{x) = Y,ei{x){V-^).. (17) 

i=l 

and is independent of A^,. The special spatial points 
selected from a discrete set of points along x, 
as well as the basis can be found using Alg. (1) which is 
described in the next section. 

We note that the kernels IC{Xx] x) appear explicitly on 
the right hand side of (16). Because of this, we have to be 
able to directly evaluate the kernel at the empirical inter¬ 
polation nodes Fortunately this is possible in 

this case as we have closed form expressions for the ker¬ 
nels. If the kernels were solutions to ordinary or partial 
differential equations that needed to be evaluated numer¬ 
ically then using the empirical interpolant becomes more 
challenging, however this is not required here (see, e.g., 
[30, 32, 33] for applications of the empirical interpola¬ 
tion method to ordinary and partial differential equation 
solvers). 

C. The Empirical Interpolation Method Algorithm 
(Offline) 

The empirical interpolation method algorithm solves 
(16) for arbitrary A^,. While it would be possible in prin¬ 
ciple to use arbitrary basis functions, such as Lagrange 
polynomials which are common in interpolation problems 
[34, 35], we take a different approach that uses only the 
information contained in the kernels themselves. We will 
take as our basis a set of M judiciously chosen kernels 
sampled at points on the parameter space where 

M is equal to the number of basis elements in (16). Be¬ 
cause the kernels vary smoothly with A^, a linear combi¬ 
nation of the basis elements will give a good approxima¬ 
tion to IC{Xx]x) for any parameter value [30]. We can 
then build an interpolant using this basis by matching 
JC{Xx;x) to the span of the basis at a set of M interpola¬ 
tion nodes The empirical interpolation method 

algorithm, shown in Alg. (1), provides both the basis and 
the nodes. 

The empirical interpolation method algorithm uses a 
greedy procedure to select the reduced basis elements 
and interpolation nodes. With the greedy algorithm, the 
basis and interpolant are constructed iterative whereby 
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the interpolant on each iteration is optimized according 
to an appropriate error measure. This guarantees that 
the error of the interpolant is on average decreasing and 
- as we show in section III D - that the interpolation error 
decreases exponentially quickly. We follow Algorithm 3.1 
of [36] which is reproduced in Alg. (1). 


The first input to the algorithm is a training space (TS) 
of kernels - distributed on the parameter space - and 
the associated set of parameters. This training space is 
denoted by T = {A^ , ^)}iLi should be densely 

populated enough to represent the full space of kernels as 
faithfully as possible. Hence it is important that 1 ^ A^. 
The second input is the desired maximum error of the 
interpolant e. We find that the norm is a robust 
error measure for the empirical interpolant and hence e 
corresponds to the largest tolerable difference between 
the empirical interpolant and any kernel in the training 
set T. 


The algorithm is initialized on steps 3 and 4 by set¬ 
ting the zeroth order interpolant to be zero, and defining 
the zeroth order interpolation error to be infinite. The 
greedy algorithm proceeds as follows: We identify the 
basis element on iteration i to be the IC{Xx]x) G T that 
maximizes the norm with the interpolant from the 
previous iteration, X^_i[/C](Aa,; x). This is performed in 
steps 7 and 8. On step 9 we select W, the i^^ interpola¬ 
tion node, by selecting the position at which the largest 
error occurs, and adding that position to the set of inter¬ 
polation nodes. By definition, the interpolant is equal to 
the underlying function at the interpolation nodes and 
so the error at Xi - which is the largest error on the 
current iteration - is removed. On step 10 we normalize 
the basis function. This ensures that the matrix (14) is 
well conditioned. On steps 11 and 12 we compute (14) 
and (17), which are used to construct the empirical in¬ 
terpolant (16). Finally, on step 13 we compute the inter¬ 
polation error ai between the interpolant on the current 
iteration X^[/C](Aa,; x) and lC{Xx] x) G T as in step 7. The 
procedure is repeated until < e. 


Once the interpolant for JC{Xx;x) is found, the equiv¬ 
alent interpolant for IC{Xy]y) is obtained trivially from 
Xm[IC]{Xx;x) by setting x ^ y. 


Algorithm 1 Empirical Interpolation Method Algorithm: 
The empirical interpolation method algorithm builds an in¬ 
terpolant for the kernels (6) iteratively using a greedy proce¬ 
dure. On each iteration the current interpolant is validated 
against a “training set” T of kernels and the worst interpo¬ 
lation error is identified. The interpolant is then updated so 
that it describes the worst-error point perfectly. This is re¬ 
peated until the worst error is less than or equal to a user 
specified tolerance e. 

1: Input: T = {A^ and e 

2: Set z = 0 
3: Set Xo[/C](Ax; a:) = 0 
4: Set (7o = oc 
5: while ai > e do 
6: i^i-\-l 

7: A^ = arg max ||/C(Aa,; a:) - li-i[JC]{Xx; x)\\loo 


8: 

6 (a;) = 

--tCiXUx) 

9: 

Xi = argmax ^i(x) - li-i[^i]{x)\ 

X 

10: 

ei(x) = 

- u(Xi)-Xi_ilu](Xi) 

11: 

Vlm = 

ei{Xm) 1 <i,m <i 

12: 

Bm (^) 

= ei(x) l<hrn<i 

13: 

ai = max /C(Aa;;a;) — li[JC]{Xx; x)\\l^ 

Aj; G7^ 

14: 

end while 

15: 

Output: 

Interpolation matrix {Bj(x)}jL 


polation nodes {Xj}^i. The equivalent interpolant 
for JC{Xy;y) is obtained trivially from {Bj {x)}jLi and 
{Xj}j^i by setting x ^ y and X ^ Y. 


D. Error Bounds on the Empirical Interpolant 

Before we proceed to demonstrate the utility of the 
empirical interpolant for quickly evaluating (5) we briefly 
remark on some of the error properties of the empirical 
interpolation method. A more detailed error analysis of 
the empirical interpolant can be found in [21]. For our 
purposes the empirical interpolant possess a highly de¬ 
sirable property, namely exponential convergence to the 
desired accuracy e. It can be shown [23, 36] (though we 
do not do so here) that there exists constants c > 0 and 
a > log (4) such that for any function / the empirical 
interpolant satisfies 

||/-2:m[/]||loo (18) 

This states that under the reasonable assumption that 
there exists an order M interpolant that allows for ex¬ 
ponential convergence, then the empirical interpolation 
method will ensure that we converge to this interpolant 
exponentially quickly. This is an important property as 
it means that the order of the interpolant, M, tends to 
be small for practical purposes. In addition, because the 
quantity on the right hand side is set to 
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a user specified tolerance e then we can set an a priori 
upper bound on the worst-fit of the interpolant. How¬ 
ever, one must still verify that the interpolant describes 
functions outside the training a postiori, though the error 
bound should still be satisfied provided that the training 
set was dense enough. In fact, it can be shown [23] that 
the empirical interpolation method is a near optimal so¬ 
lution to the Kolmogorov n-width problem in which one 
seeks to find the best M-dimensional (linear) approxima¬ 
tion to a space of functions. 

It is important to recall that in this paper we are in¬ 
terpolating the integral kernels (6) which are a function 
of six free parameters A^,: two indices n and n' and two 
complex beam parameters Px and q'^. Had we not used 
the EIM, we would have had to find an alternative way 
of expressing the Aa^-dependent coefficients in (12). Con¬ 
sider, for example, a case in which we had used tensor- 
product splines to describe the coefficients: Using a grid 
of just ten points in each of the six parameters in A^^ 
would result in an order 10^ spline which would surely 
be computationally expensive to evaluate. Furthermore, 
there would be no guarantee of its accuracy or conver¬ 
gence to a desired accuracy. 

E. Reduced order quadrature (Online) 

Substituting the empirical interpolant (16) into (9) 
gives the ROQ, 

^nm,n'm'i^Qx') Qx") Qy! Qy) 

M M 

^ ^ Wki /C(A,; Xfe) IC{\y; Yi) , (19) 

k=l 1=1 

with the ROQ weights uj^i given by: 

L L 

^kl = F '^WijA{xi,yj)Bk{xi)Bi{yj). (20) 

i=i j=i 

The ROQ form of the coupling coefficient enables fast 
online evaluations of the coupling coefficients. Note that 
because only operations are required to perform the 
double sum (19) we expect that the ROQ is faster than 
the traditional L^-term Newton-Cotes integration by a 
factor of provided that M < L. We expect in 

practice that M L due to the exponential convergence 
of the empirical interpolation method. 

The number of operations in (19) can be compressed 
further still due to the separability of the empirical in¬ 
terpolant (16) into beam parameters A^, and spatial po¬ 
sition X that allows us to exploit the spatial symme¬ 
try in the HG modes. The HG modes exhibit spatial 
symmetry/antisymmetry under reflection about the ori¬ 
gin. Hence it is useful to split the x and y dimensions 
into four equally sized quadrants and perform the ROQ 
in each quadrant separately. For example, when a HG 
mode is symmetric between two or four of the quadrants 


then only two or one set(s) of coefficients {/C(Aa,; 
needs to be computed (and likewise for . 

This will speed up the computation of the ROQ (19) by 
up to a factor of four. Hence, in practice we need only 
build the El over one half plane for either positive or 
negative values of x (or equivalently y)] we derive the 
basis spanning the second half-plane by reflecting the 
basis about the origin. To ensure that this symmetry 
is exploitable the data points of the map must be dis¬ 
tributed equally and symmetrically about the beam axis 
{{x^y) = (0,0)). Those points that lie on the x and y 
axes must also be weighted to take into account they 
contribute to multiple quadrants when the final sum is 
computed. In the cases where the map data points are 
not correctly aligned we found that bilinear interpolation 
of the data to retrieve symmetric points did not introduce 
any significant errors. However, higher-order interpola¬ 
tion methods can introduce artefacts to the map data. 

IV. EXEMPLARY CASE: NEAR-UNSTABLE 
CAVITIES AND CONTROL SIGNALS 

There are several scenarios when modelling tools can 
benefit heavily from the ROQ method, of particular in¬ 
terest are cases where the simulation time is dominated 
by the integration time of the mirror surface maps. One 
such example is an investigation into the feasibility of 
upgrading the LIGO interferometers with near-unstable 
arm cavities. The stability of a Fabry-Perot cavity is de¬ 
termined by its length L and radius of curvature (RoG) 
of each of its mirrors and can be described using the pa¬ 
rameter:: 

g = {1 — L / Rcii^)(l — L / Rcqi^). (21) 

with 0 < ^ < 1 defining the stable region. Near-unstable 
cavities are of interest because they result in larger beam 
sizes on the cavity mirrors (see also figure 3) which re¬ 
duces the coating thermal noise [11], one of the limiting 
noise sources of the detector. One negative aspect of such 
near-unstable cavities is that the transverse optical mode 
separation frequency approaches zero as ^ ^ 0 or 1. The 
mode separation frequency determines the difference in 
resonance frequency of higher-order modes with respect 
to the fundamental mode. Thus with a lower separation 
frequency any defect in the cavity causing scattering into 
HOMs is suppressed less and can contaminate control 
signals for that cavity and couple extra noise into the 
GW detectors output The optimal cavity design must 
be determined as a trade-off between these degrading ef¬ 
fects and the reduction in coating thermal noise. This is 
a typical task where a numerical model can be employed 


^ Another potential problem is additional clipping or scattering of 
the beam on the mirrors due to the larger beam sizes which can 
result in increased roundtrip losses of the arm cavity. 
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(a) Cavity scan as ITM and ETM RoC varied 
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(b) The relative error between ROQ and Newton-Cotes. 


Figure 2: Modeled LIGO cavity scan as the RoC of the ITM and ETM are varied to make the cavity increasingly more 
unstable. This simulation was run for Omax = 10 and includes clipping from the hnite size of the mirrors and surface 
imperfections from the ETM08 and ITM04 maps. Eigure 2a shows how the amount of power scattering into HOM changes as 
g ^ 1. Also visible here is the reduction in the mode separation frequency with increasing instability. The contribution of the 
TEMoo mode has been removed to make the HOM content more visible. The reduced basis was built for mode order O — 14, 
to reduce errors, see hgure 7. The difference in this result when using ROQ compared to Newton-Cotes is shown in 2b. 
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Figure 3: The beam size on the ITM and ETM of a LIGO 
cavity as a function of cavity stability parameter as the 
mirror RoCs are tuned. 


to search the parameter space. In this case each point 
in that parameter space corresponds to a different beam 
size in the cavity which forces a re-computation of the 
scattering matrices on the mirrors. Thus the new algo¬ 
rithm described in this paper should yield a significant 
reduction in computing time. 

In this section we briefiy summarise the results from 
the simulations and in the following section we provide 
the details of setting up the model and give an analysis 


of the performance of the ROQ algorithm. We have im¬ 
plemented the ROQ integration in our open-source sim¬ 
ulation tool Finesse and use the official input parame¬ 
ter files for the LIGO detectors [20]. Below we show the 
preliminary investigation of the behaviour of a single Ad¬ 
vanced LIGO like arm cavity with a finesse of 450, where 
the mirror maps for the mirrors ETM08 and ITM04 ^ 
were applied to the high reflective (HR) surfaces. Note 
that we do not report the scientific results of the sim¬ 
ulation task which will be published elsewhere. This 
example is representative for a class of modelling per¬ 
formed regularly for the LIGO commissioning and de¬ 
sign and provides us with a concrete and quantitive setup 
to demonstrate the required steps to use the ROQ algo¬ 
rithm. 

Modelling the LIGO cavity for differing stabilities in¬ 
volves varying the RoG of both the ITM and ETM. 
The resulting change in w{z) at each surface means the 
scattering matrices will need to be recomputed for each 
state we choose. To view the HOM content in the cav¬ 
ity created by the scattering a cavity scan can be per¬ 
formed, displacing one of the cavity mirrors along the 


^ The nominal radius of curvatures of ETM08 and ITM04 are 
1934m and 2245m respectively. The optical properties of these 
mirrors were taken from [25]. 
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Figure 4: The Pound-Dr ever-Hall error signal for the 
LIGO cavity modelled in figure 2. A significant change in 
zero-crossing position and shape can be seen as the stability 
of the cavity is reduced {g 1). 


cavity axis on the order of the wavelength of the laser 
light, A = 1064 nm, to change the resonance condition 
of the cavity. We have performed the simulations using 
O = 10 with Newton-Cotes integration and our ROQ 
method. The results for cavity scans at different RoCs 
are shown in figure 2a. The dominant mode is the funda¬ 
mental TEMoo whose resonance defines the zero tuning, 
the power in the TEMqo mode has been removed from 
this plot to better show the lower power HOM content. 
Eor more stable cavities (at the bottom of the plot in fig¬ 
ure 2a) the HOMs are well separated and not resonant at 
the same time as the TEMqo- As the RoC is increased, 
the stability is reduced and the HOMs can be seen to 
converge and eventually become resonant at a tuning of 
0. At a stability of g ~ 0.98 the cavity mode begins to 
break down significantly and many modes become reso¬ 
nant. The effect of this on a sensing and control signal 
used for a Pound-Drever-Hall control system is shown in 
figure 4, where for increasingly unstable cavities the error 
signal becomes degraded, showing an offset to the nomi¬ 
nal zero crossing, a reduced slope and overall asymmetry 
around the centre. The complete investigation into the 
feasibility of such cavities is beyond the scope of this 
paper, it includes amongst other issues the quantitative 
comparison of the control noise from the degradation of 
the control signals with the reduced thermal noise. The 
simulation task described above is sufficient to provide a 
test case for our ROQ method. 


V. APPLICATION AND PERFORMANCE OF 
NEW INTEGRATION METHOD 

In this section we provide a detailed and complete 
recipe for setting up and using the ROQ for the LIGO 
example, using Finesse and Pykat, and discuss the per¬ 
formance, in terms of speed and accuracy, of our method. 


The description should be sufficient for the reader to im¬ 
plement our method for their own optical setup. In this 
section we include the costly “offline” procedure of build¬ 
ing the empirical interpolant for completeness. As part 
of the ongoing code development of Finesse and Pykat 
we intend to pre-generate the empirical interpolants, suit¬ 
able for a wide range of problems, so that the typical user 
should not need to run the costly offline building of the 
interpolant into their simulations. 


A. Computing the ITM and ETM Empirical 
Interpolants 


Firstly the range of beam parameters for the simula¬ 
tion must be determined. Once this is known a training 
set can be constructed and the empirical interpolant can 
be computed. The surface distortions that are of inter¬ 
est are those on the HR surfaces of a LIGO arm cavity 
mirror. We will require two FIs, one for the ITM HR sur¬ 
face and one for the ETM HR surface due to the differing 
beam parameters at each mirror. The beam parameter 
range that the training sets should span are determined 
by varying the radius of curvature of the ITM and ETM 
to include the range of cavity stabilities which we want 
to model. The beam parameter ranges are shown in fig¬ 
ure 5. The required ranges for the ETM are 4.7m < wq < 
12.0mm and 2.11km < z < 2.20km and for the ITM 
4.7mm < wq < 12.0m and —1.88km < 2 ; < —1.79km, up 
to a maximum optical mode order of O = 20, Netwon- 
Cotes degree of 6, L = 1199. For this example we fix the 
maximum tolerable error of the empirical interpolant to 
e = 10-^1 

Using these ranges the method described in sec¬ 
tion III A can be used to produce the FIs. The offline 
computation of the basis can have significant computa¬ 
tional cost. For very wide parameter ranges the memory 
required to store the training sets can quickly exceed that 
of typical machines. For the above parameters, with 100 
sample points each in the wq and z range, up to O = 14 
and e = 10“^^ approximately 7GB of memory was re¬ 
quired. Running this method on machines with less mem¬ 
ory is possible by storing the training set on a hard drive 
using a suitable data storage format such as HDF5 for 
access. Computation time of the empirical interpolant is 
then limited by the read and write times of the media. 
Using a MacBook pro 2012 model which contains a 2.7 
GHz Intel core i7 with 8GB of RAM generating the ITM 
and ETM reduced basis and empirical interpolant takes 
~ 4 hours each. The number of elements in the final 
reduced basis for the ITM and ETM were N = 30 and 
N = 29 respectively. In figure 8 the convergence of the 
empirical interpolant error with respect to the acceptable 
empirical interpolant error. One can see that the El error 
converges exponentially as described in section HID. 
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Figure 5: Range of beam parameters needed to model a change in curvature from Om to 90 m at the ITM and the ETM. In 
order to utilise the ROQ to cover this parameter space, the empirical interpolant needs to be constructed using a training set 

made from kernels (6) densely covering this space. 


B. Producing the ROQ weights 

Once the empirical interpolant has been computed for 
both ITM and ETM HR surfaces the ROQ weights (20) 
can be computed by convoluting the mirror maps with 
the interpolant. The surface maps that we have chosen 
are the measured surface distortions of the (uncoated) 
test masses currently installed at the LIGO Livingston 
observatory, shown in figure 1. The maps contain L « 
1200 samples and we can expect a theoretical speed-up 
of ^ 1200^/30^ = 1600 from using ROQ over 

Newton-Cotes. These maps include an aperture, M, and 
the variation in surface height in meters, z{x,y). Thus 
to calculate the HOM scattering on reflection from one 
of these mirrors with (5) the distortion term is: 

A{x,y)=A{x,y)e^'^’^^^-<y^ (22) 

where A{x, y) is 1 if ^/x^ + y^ < 0.16in and 0 otherwise, 
and k is the wavenumber of the incident optical field. 

Using (22) with equation (20) (with a Newton-Cotes 
rule of the same degree the empirical interpolant was gen¬ 
erated with) the ROQ weights can be computed for each 
map shown in figure 1. This computational cost is pro¬ 
portional to the number of elements in the El, M, and the 
number of samples in the map, L^. For the LIGO maps 
this takes ~ 10s on our 2012 MacBook Pro. The result¬ 
ing ROQ rule for the maps can be visualised as shown in 
figure 6: the amplitude of the ROQ weights map out the 
aperture and the phase of the weights varies for different 
maps because of the different surface structure. The com¬ 
putation of these ROQ weights need only be performed 
once for each map, unless the range of beam parameters 
required for the empirical interpolant are changed. 

We verify that the process of generating the ROQ rule 
has worked correctly by computing the scattering ma¬ 
trices with ROQ and Newton-Cotes across the parame¬ 
ter space. We compute k{q; ETM07) with Omax = 10 
using ROQ and then again using Newton-Cotes integra¬ 
tion. Computing the relative error between each element 
of these two matrices the maximum error can be taken 
for q values spanning the requested q parameter range. 


Figure 9 shows how the final error of the El, ctm, propa¬ 
gates into an error in the scattering matrix. This shows 
the maximum (solid line) and minimum (dashed line) er¬ 
rors for any element in the scattering matrix between 
the two methods. From this it can be seen that building 
a more accurate empirical interpolant results in smaller 
maximum errors in the scattering matrix. Now, using 
the most accurate reduced basis the maximum relative 
error is shown in figure 7 over the q space, where the 
white dashed box shows the boundaries of the parame¬ 
ters in the training set. Overall the method successfully 
computed a ROQ rule that accurately reproduced the 
Newton-Cotes results for scattering up to O = 10. In 
should be notes that the largest errors, e.g. as seen in 
figure 9, do not represent the full parameter space but 
occur only at smallest 2 : and largest wq. It was also found 
that building a basis including a higher maximum HOM, 
for example basis of order 14 for scattering computations 
up to order 10, significantly improved the accuracy of 
the ROQ. Using an reduced basis constructed for order 
14 rather than order 10 only increased the number of el¬ 
ements in the basis by 2, thus not significantly degrading 
any speed improvements. It can also be seen in figure 7 
that ROQ extrapolates beyond the originally requested 
q parameter space and does not instantly fail for evalu¬ 
ations outside of it. A gradual decrease in the accuracy 
can be seen when using larger wq values. 


C. Performance 

The time taken to run these Finesse simulations as 
O is increased is shown in figure 10 demonstrating how 
much more efficient it is to use ROQ over Newton-Cotes 
for the computation of scattering matrices. We also show 
for reference the computation time when no scattering 
from surface maps is included to give the base time it 
takes to run the rest of the Finesse simulation. The 
overall speed-up achieved can be seen in figure 11, reach¬ 
ing 2700 times faster to run the entire simulation 
at O = 10. The overall speed-up then begins to drop 
slightly as the base time taken to run the rest of Finesse 
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Figure 6: Absolute and argument values of the ROQ weights (20) generated for each of the maps as shown in figure 1. Here 
the final quadrature rule can be visualized. The top plots show the absolute value: The size of the point is proportional to |u;| 
and the center of each point lies on a specific empirical interpolation node in the x-y plane (Xi,Yj) (c.f. (13) and (19)). The 
bottom plots show log^Q{arg{w)). The dashed line on each plot shows the mirror surface boundary; outside the boundary the 
mirror maps are equal to zero. We note that there are non-zero ROQ weights associated with points in the region where the 
mirror maps are zero. While this may be counter intuitive, it is a consequence of the fact that the empirical interpolant nodes 
lie within the full x-y plane and, that they are constructed without any knowledge of the mirror maps: the weights still 
receive no contribution from the region where A{x,y) = 0 as this region does not contribute to the sum in (20). However, the 
ROQ uses information about the kernels (6) over the entire region, including where A(x,y) = 0. 
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Figure 7: Maximum relative error between the scattering matrices computed for the ETM07 surface map, with ROQ and 
using Newton-Cotes, for mode orders up to Omax = 18. The dashed white area represents the beam-parameter region over 
which training sets were generated. The subplots illustrates how using an ROQ built for a larger Omax scattering reduces the 
maximum error significantly. Also shown is that the ROQ is valid over a larger parameter range than what it was initially 
generated for, implying that the empirical interpolant can be used for extrapolation in a limited parameter region outside the 

initial range indicated by the white dashed box. 


becomes larger. The dashed line in figure 11 shows the 
speed-up if this base time is removed, again showing an 
impressive speed-up peaking at 4000 times faster. 

Using ROQ enables us to perform such modelling tasks 
with a far greater efficiency. Running the model to com¬ 
pute the output seen in figure 2a required computing 100 
different scattering matrices for the various changes in 


RoC. This took 20.5 hours to compute with Newton- 
Cotes and 18 minutes with ROQ The difference in the 


^ Note that the effective speed-up in this case is less than the 
large values in figure 11 because here we have included the total 
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Basis size 

Figure 8: El error as a function of the number of basis 
elements selected by the greedy algorithm (Alg. 1) for the 
example described in section V A. As expected from the 
error analysis in section HID, the empirical interpolant error 
displays exponential convergence with the basis size. 



El tolerance, e 

Figure 9: Relative error in the scattering matrices 
computed using the ROQ and Newton-Cotes integration 
(with Omax — 10) as a function of the empirical interpolant 
tolerance e. The empirical interpolant was built for 
maximum coupling Omax = 14. The error is the minimum 
(dashed lines) and maximum (solid lines) over the parameter 
space with which the empirical interpolant was built for, 
thus represents the worst and best case scenarios. The 
largest errors are independent of the map data and occur on 
couplings coefficients which couple the higher order modes 
included in the empirical interpolant. 


final result between ROQ and Newton-Cotes is shown in 
terms of relative error in figure 2b. We have prepared the 
ROQ input for this example such that the error is signif¬ 
icantly lower than Ippm (relative error of 10“^) thereby 
showing that ROQ can be both much faster and still suf¬ 
ficiently accurate. 


runtime of the simulation. This includes the initialisation and 
running of the other aspects of Finesse which took 17 minutes. 
The actual time taken for just the ROQ calculation is 30 s thus 
a speed-up in the ROQ vs Newton-Cotes is 2500. 
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Figure 10: Time taken to run Finesse to model the steady 
state optical fields in an LIGO cavity with surface maps on 
both the ITM and ETM HR surfaces. The timing of running 
the entirety of Finesse is used—rather than just the core 
method—because there are additional speed improvements 
from having to read and handle significantly less data points, 
from the L x L maps down to M x M ROQ weights. Smaller 
data fits into processor caches better and also reduces disk 
read times. This plot compares a single computation of the 
scattering matrices with ROQ (19) and Newton-Cotes. The 
case with no maps used is also shown to illustrate how much 
time is spent in Finesse doing calculations not involving 
maps, which now becomes the dominant computational cost 
when using ROQ. A significant improvement is also found 
for order zero where only one scattering integral need be 
calculated; this is partly time saved from having to read 
larger data from the disk and manipulating it in memory. 
The preprocessing is unavoidable as the Finesse can except 
many different types of map, thus it cannot be optimised at 
runtime until it know what it is dealing with. ROQ helps 
here as it removes this preprocessing step so it need only 
happen once. The ROQ preprocessing happens during the 
computing of the ROQ weights (20). This is a one time cost 
for each map for a particular El donqe outside of Finesse, 
thus isn’t included in this timing. The computational cost of 
this is on the order of 5s for each map for the reduced basis 
used in this example. 


VI. CONCLUSION 

Numerical modelling of optical systems plays a vital 
role for the design and commissioning of precision inter¬ 
ferometers. The typical use of the simulation software 
in this area requires rapid iterations of many simula¬ 
tion runs and manual fine tuning as modelling progresses, 
which is not well suited for large computer clusters. The 
scope of current investigations is often limited by the re¬ 
quired computation time and thus the development of 
fast and flexible tools is a priority. Current problems 
in precision interferometers, such as LIGO, involve the 
investigation of laser beam shape distortions and their 
effect on the interferometer signal. Frequency-domain 
simulations using Gaussian modes to describe the beam 
properties have emerged as fast and flexible tools. How¬ 
ever, the computation of the scattering matrix for mir- 
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Figure 11: The speed-up achieved using ROQ compared to 
Netwon-Cotes as a function of mode order using the timing 
values in figure 10. The dashed line shows the speed-up if 
the time for initialisation and post-processing is subtracted 
from both times for Newton-Cotes and ROQ. This 
demonstrates the improvements just for the computational 
cost relating to map scattering calculations. Simulations 
that have a larger computational cost relating to features 
not related to scattering will show a smaller speed-up. For 
example, the simulations results shown in figure 2 have a 
total simulation speed-up of ~ 80 but the scattering 
calculation was reduced from 20.5 Hours ^ 30 s. 


ror surface distortions—effectively an overlap integral of 
measured surface data with Hermite-Gauss modes—has 
shown to be a limiting factor in improving the compu¬ 
tational speed of such tools. A significant reduction in 
computational time of current numerical tools is required 
for more efficient in-depth modelling of interferometers 
including more realistic features such as clipping, optical 
defects, thermal distortions and parametric instabilities. 

In this work we have demonstrated how the empiri¬ 
cal interpolation method can be used to generate an op¬ 
timised quadrature rule for paraxial optical scattering 


calculations, known as a reduced order quadrature. Our 
method removes the prohibitive computational cost of 
computing the scattering by speeding up the calculation 
of the steady state optical fields in a LIGO arm cavity by 
up to a factor of 2750 times, reducing simulation times 
from days to minutes. Using an exemplary simulation 
task of near-unstable arm cavities for the LIGO interfer¬ 
ometers we have demonstrated that our method is both 
accurate and fast for a typical modelling scenario where 
imperfections in the interferometer have a significant im¬ 
pact on optical performance. We have provided a com¬ 
plete recipe to recreate and use the new algorithm and 
provide an open source implementation in our general- 
purpose simulation software Finesse. Importantly, the 
reduced order quadrature integration method is generic 
and can be applied to any optical scattering problem for 
any surface distortion data. 
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